/*************************************************************************
Purpose: This program analyzes timeseries in universal plot format.
         It differs from grace_input_timeseries.cpp because that program
         uses input data as input, while this one reads universal 2D plots.
**************************************************************************
Input:   Reads output_filenames.txt in output/ to find files with 
         timeseries.
**************************************************************************
Output:  2D plot files
**************************************************************************
Copyright 2010 Bryan Killett

This file is part of GAIA.

GAIA is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.

GAIA is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.

You should have received a copy of the GNU General Public License
along with GAIA. If not, see <http://www.gnu.org/licenses/>.
*************************************************************************/
//Header file "definitions.hpp" declares many functions, constants, datatypes.
//In order to use the functions declared in definitions.hpp, must link
//with object code that includes functions.cpp.
#include "definitions.hpp"

typedef struct { 

  vector <double> real,imag;
} complex_s;

int custom_gsl_fft_halfcomplex_unpack (const double halfcomplex_coefficient[],
                                       complex_s &transform,
                                       const size_t n);

int custom_gsl_fft_halfcomplex_unpack (const double halfcomplex_coefficient[],
                                       complex_s &transform,
                                       const size_t n){

  size_t i;//first time I've used this type!
  
  //Just in case.
  transform.real.clear(); transform.imag.clear();
  transform.real.resize(n, 0.0); transform.imag.resize(n, 0.0);
  
  transform.real.at(0) = halfcomplex_coefficient[0];
  transform.imag.at(0) = 0.0;
  
  for(i = 1; i < n - i; i++){
    double hc_real = halfcomplex_coefficient[(2 * i - 1)];
    double hc_imag = halfcomplex_coefficient[(2 * i)];
    transform.real.at(i) = hc_real;
    transform.imag.at(i) = hc_imag;
    transform.real.at(n-i) = hc_real;
    transform.imag.at(n-i) = -hc_imag;
  }
  
  if(i == n - i){
    transform.real.at(i) = halfcomplex_coefficient[(n - 1)];
    transform.imag.at(i) = 0.0;
  }
  
  return 0;
}

//This function is defined at the bottom of this file.
int timeseries_analysis(results_s &input,
                        results_s &output);

int main(int argc, char *argv[]){

  string basefolder(gaiafolder());
  string file;
  results_s input, output;
  plot_options_s input_plot_options, output_plot_options;
  input_s dummy_input;
  grid_s dummy_grid;
  ancillary_s dummy_ancillary;
  //These variables are used to display run time of program after execution.
  time_t start_time,end_time;
  timespan_conversion_s timespan_conversion;
  string base;//This string forms the base part of the output filename.
  
  start_time = time(NULL);//Record time at start of program.

  //Use current time and command line args to create output folder.
  create_output_folder(argc,argv,output_plot_options);

  //Each program should manually set this value.
  //See plot_options in definitions.hpp for more details.
  input.options.output_choice  = 2;
  output.options.output_choice = 2;
  //Define input folder.
  input_plot_options.inputfolder = "output/maps/1249576470 - autoFES k1 and its load file ONLY not zeroed vs. masconFES k1 only, custom grid 60N, 330E, TIME SERIES/time series/0 offset, 20 years/";

  //   1 - Calculate FFT of input data, return as real and imag.
  //   2 - Calculate FFT of input data, return as amp and phase.
  output.options.type=2;
  
  //Control plotting FFTs (in hours).
  output.options.min_period = 22;
  output.options.max_period = 25;
  
  //Load files in the input folder.
  load_output(input, input_plot_options, dummy_grid);

  if(input.options.output_choice != 2) cout<<"WARNING!!! input.output_choice (should be 2) == "<<input.options.output_choice<<endl;

  timeseries_analysis(input, output);

  //Set the output filename base.
  base = "plot_data";
  
  //Now that the analysis has been performed, it's time to output the data
  //in a format that the IDL program grace_input_map.pro can turn into a 
  //plot.  Last variable is 1 because this is the final write_output.
  write_output(base, output, dummy_grid, output_plot_options,1);

  end_time = time(NULL);//Record time at end of program.
  timespan_conversion.seconds = (long)(end_time - start_time);
  cout<<"analyze_timeseries took"<<sec2human(timespan_conversion)<<"."<<endl;

  return 0;
}

int timeseries_analysis(results_s &input,
                        results_s &output){
  /**********************************************************************
  Purpose:
            This function 
  ***********************************************************************
  Input:  (variable types ending in "_s" are defined in definitions.hpp)
  ***********************************************************************
  Output: (variable types ending in "_s" are defined in definitions.hpp)
  ***********************************************************************
  Written by Bryan Killett
  First written:  2009-08-06
  Last modified:  2009-08-06
  **********************************************************************/
  long i,j,k;

  string temp_string;
  
  vector<double> temp_v_double;//Used to push_back vectors of dbl vectors
  vector<string> temp_str_double;//Used to push_back vectors of str vectors
  vector< vector<double> > temp_vv_double;//push_back vs of vs of dbl vectors!

  //Perform different types of timeseries analyses based on "type".
  switch(output.options.type){
    //////////////////////////////////////////////////////////////////////
    //////////////////////////////////////////////////////////////////////
    //1 - Calculate FFT of input data, return as real and imag.
    //2 - Calculate FFT of input data, return as amp and phase.
    case 1:
    case 2:{

      gsl_fft_real_wavetable * real;
      gsl_fft_real_workspace * work;
      
      long outi=0;//Like cntr "i" but for output (2x size of input!)

      //Loop through different plots.
      for(i=0;i<(long)input.xy.x_values.size();i++){

        //Save misc info for this REAL/AMP plot.
        temp_string = input.xy.titles.at(i);
        temp_string.append(" FFT");
        output.xy.titles.push_back(temp_string);
        output.xy.x_units.push_back("Period (hours)");
  
        if(output.options.type==1){
          output.xy.y_units.push_back("Real part of FFT (cm)");
        }
        else if(output.options.type==2){
          output.xy.y_units.push_back("Amplitude of FFT (cm/sqrt(Hz))");
        }       
        
        output.xy.legends.push_back(temp_str_double);

        //This is a brand new plot, so push_back using the new vv.
        output.xy.x_values.push_back(temp_vv_double);
        output.xy.y_values.push_back(temp_vv_double);
        
        //Save misc info for this IMAG/PHASE plot.
        temp_string = input.xy.titles.at(i);
        temp_string.append(" FFT");
        output.xy.titles.push_back(temp_string);
        output.xy.x_units.push_back("Period (hours)");
        
        if(output.options.type==1){
          output.xy.y_units.push_back("Imag part of FFT (cm)");
        }
        else if(output.options.type==2){
          output.xy.y_units.push_back("Phase of FFT (degrees)");
        }
        
        output.xy.legends.push_back(temp_str_double);

        //This is a brand new plot, so push_back using the new vv.
        output.xy.x_values.push_back(temp_vv_double);
        output.xy.y_values.push_back(temp_vv_double);
        
        //Loop through different timeseries in each plot.
        for(j=0;j<(long)input.xy.x_values.at(i).size();j++){
        
          output.xy.legends.at(outi).push_back(input.xy.legends.at(i).at(j));//REAL/AMP
          output.xy.legends.at(outi+1).push_back(input.xy.legends.at(i).at(j));//IMAG/PHASE

          long numpoints = (long)input.xy.x_values.at(i).at(j).size();
          gsl_vector * values = gsl_vector_calloc(numpoints);
          complex_s transform;

          //delta = time resolution of this particular time series (assumed uniform).
          //AND ASSUMED TO BE IN HOURS!!!
          double delta = input.xy.x_values.at(i).at(j).at(1) - input.xy.x_values.at(i).at(j).at(0);
          //New timeseries for this REAL/AMP plot.
          output.xy.x_values.at(outi).push_back(temp_v_double);
          output.xy.y_values.at(outi).push_back(temp_v_double);
          //New timeseries for this IMAG/PHASE plot.
          output.xy.x_values.at(outi+1).push_back(temp_v_double);
          output.xy.y_values.at(outi+1).push_back(temp_v_double);
          //Set up GSL temp variables.  Needs to be allocated in loop b/c
          //the size is potentially different each time.
          real = gsl_fft_real_wavetable_alloc (numpoints);
          work = gsl_fft_real_workspace_alloc (numpoints);
          for(k=0;k<numpoints;k++) gsl_vector_set(values,k,input.xy.y_values.at(i).at(j).at(k));

          //Perform an FFT on the input accels, which are real numbers.
          gsl_fft_real_transform(values->data, 1, numpoints, real, work);
          //Free up memory taken up by GSL temp vars.
          gsl_fft_real_wavetable_free (real);
          gsl_fft_real_workspace_free (work);
          custom_gsl_fft_halfcomplex_unpack(values->data, transform, numpoints);
          gsl_vector_free(values);

          //Save individual XY values for this plot.
          //Only loop to the middle of the number of input points in the
          //swath, because after that point the negative frequencies start
          //and they're just complex conjugates of the previous ones.
          //Skipping index 0 because that's the zeroth frequency component.
          double base_frequency = 1/(delta*numpoints);
          for(k=1;k<(numpoints)/2;k++){
            double current_period = 1/(k*base_frequency);
            if(current_period >= output.options.min_period and current_period <= output.options.max_period){
              if(output.options.type==1){
                //Real
                output.xy.x_values.at(outi).at(j).push_back(current_period);
                output.xy.y_values.at(outi).at(j).push_back(transform.real.at(k));
                //Imag
                output.xy.x_values.at(outi+1).at(j).push_back(current_period);
                output.xy.y_values.at(outi+1).at(j).push_back(transform.imag.at(k));
              }
              else{
                //AMP
                output.xy.x_values.at(outi).at(j).push_back(current_period);
                output.xy.y_values.at(outi).at(j).push_back(sqrt(pow(transform.real.at(k),2)+ pow(transform.imag.at(k),2)));
                
                //PHASE
                output.xy.x_values.at(outi+1).at(j).push_back(current_period);
                output.xy.y_values.at(outi+1).at(j).push_back(rad2deg*atan2(transform.imag.at(k),transform.real.at(k)));
              }
            }
          }
        }//End of loop through timeseries in this plot.
        
        outi+=2;//Real and imaginary plots.
      }//End of loop through different plots.
      break;
    }//End of cases 1,2.
  }//End of switch-case-break  

  return 0;
}
